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We present the extension of our COCAL - Compact Object CALculator - code to compute general-relativistic 
initial data for binary compact-star systems. In particular, we construct quasiequilibrium initial data for equal- 
mass binaries with spins that are either aligned or antialigned with the orbital angular momentum. The Isenberg- 
Wilson-Mathews formalism is adopted and the constraint equations are solved using the representation formula 
with a suitable choice of a Green’s function. We validate the new code with solutions for equal-mass binaries 
and explore its capabilities for a wide range of compactnesses, from a white dwarf binary with compactness 
~ 10~^, up to a highly relativistic neutron-star binary with compactness 0.22. We also present a comparison 
with corotating and irrotational quasiequilibrium sequences from the spectral code LORENE [Taniguchi and 
Gourgoulhon, Phys. Rev. D 66, 104019 (2002)] and with different compactness, showing that the results 
from the two codes agree to a precision of the order of 0.05%. Finally, we present equilibria for spinning 
configurations with a nuclear-physics equation of state in a piecewise polytropic representation. 


I. INTRODUCTION 

With a compactness slightly smaller than that of a black 
hole, neutron stars are most probably nature’s ultimate com¬ 
pact matter configuration before gravitational collapse and 
black-hole formation. As such, they present an invaluable tool 
to astrophysicists in order to study a plethora of problems and 
test the limits of existing knowledge, from general relativity, 
via the emission of gravitational waves, to nuclear physics, via 
the input on the equation of state of nuclear matter [2-6] . For 
example, the leading (but not unique) candidate to explain one 
of the most luminous explosions in the universe, the so-called 
short gamma-ray bursts [7, 8] (see [9] for a recent review) 
is the merger of two neutron stars (or of one neutron star and 
one black hole) with the subsequent formation of a black hole, 
an accretion torus, and a jet structure of ultrastrong magnetic 
field [10, 11]. Yet another example has to do with the produc¬ 
tion site of the heaviest elements in the universe through the 
so-called, rapid neutron capture (r-process) [12-17]. 

Central to the processes described above [18, 19] is a bi¬ 
nary neutron star system which in addition constitutes a prime 
source of gravitational waves for ground-based laser interfero¬ 
metric gravitational-wave detectors such as LIGO, Virgo, KA- 
GRA, and ET [20-24]. The advanced generation of these 
detectors will become operational in a few years, and they 
will be able to observe a volume of the universe a thousand 
times more than their predecessors. According to present es¬ 
timates [25] it may be possible to detect ^ 1 — 100 events per 
year, making the study of such systems an important step to¬ 
ward a practical verification of general relativity in the strong 
field regime, as well as an exploration of its limits. At the 


same time, gravitational-wave observations are expected to 
constrain the neutron-star equation of state [26-34]. 

The broad-brush picture for the two-body problem in gen¬ 
eral relativity can be divided into three phases: the inspiral, 
the merger, and the ring-down, with each one having its own 
methods and tools of investigation. The purpose of this work 
lies in the interface between the first and the second phases, 
the so-called quasiequilibrium stage, and the solutions pre¬ 
sented are meant as “snapshots” at particular instants of the 
binary system. The purpose is twofold: on the one hand to 
provide initial data for the simulation of the merging phase 
and, on the other hand, to provide evolutionary information 
about the system studied by constructing quasiequilibrium se¬ 
quences. In this way, we can learn how much the star shape 
is deformed as the orbit shrinks, where an instability sets in, 
where the mass-shedding limit is and what the angular veloc¬ 
ity of the system is there. All of this information can be com¬ 
puted with a modest computational infrastructure, thus allow¬ 
ing for the exploration of a wide parameter space. 

Because the orbital decay time scale due to gravitational 
wave emission is much shorter than the synchronization time 
scale due to the neutron star viscosity, it is unlikely that the 
two stars will be tidally locked before merger [35, 36]. For 
such slowly rotating configurations, the assumption of an irro¬ 
tational flow is physically reasonable and mathematically sim¬ 
ple to impose. An irrotational flow is also called a potential 
flow, since the fluid velocity is the gradient of a potential [37]. 
A formalism to compute initial data within a conformally 
flat geometry, the so-called Isenberg-Wilson-Mathews (IWM) 
formalism [38, 39], was presented in Refs. [40-43], and nu¬ 
merical implementations for a variety of physical assumptions 
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have been discussed by several groups [1, 44-53]. Nonconfor- 
mally flat formulations, have also been implemented [54-59], 
where the full system of Einstein equations is being solved. 
These more computationally expensive solutions are expected 
to respect the circularity of an orbit better than the ones com¬ 
ing from conformally flat initial data, which in addition seem 
to suppress tidal effects as the compactness of the stars in¬ 
creases. A conformally flat geometry can still be used to pro¬ 
duce low eccentricity initial data if one uses ideas similar to 
those applied to the binary black hole problem [60], as they 
were implemented in [61, 62]. 

As it was pointed in Ref. [63], the double pulsar PSR 
J0737-3039, has one of its stars reaching the merging epoch 
with a spin of ~ 27 ms, hence in a state that cannot be consid¬ 
ered irrotational. Since we only know less than a dozen binary 
systems [64], and hundreds of millisecond pulsars, it is rea¬ 
sonable to simulate arbitrary spinning binary neutron stars and 
assess the impact that the stellar spin has on the gravitational- 
wave signal. Initial data for binary neutron stars with interme¬ 
diate (and arbitrary) rotation states are more difficult to cal¬ 
culate, since there is no self-consistent scheme to incorporate 
the fluid equations with the rest of the elliptic gravitational 
equations. Various schemes have been proposed recently in 
Refs. [63, 65-68], which introduce some additional approxi¬ 
mations, while evolution of spinning binaries have been per¬ 
formed in [68-70]. 

In this work we continue the COCAL program for com¬ 
puting equilibrium configurations of single [71, 72] and bi¬ 
nary systems building on the infrastructure introduced in Refs. 
[73-75] for binary black holes. Here, we describe how to cal¬ 
culate initial data for binary stars, concentrating on neutron 
stars. The ability to compute configurations with a wide range 
of compactness was one of the goals of this work. At present, 
COCAL makes use of a piecewise polytropic description 
to represent the equation of state (EOS), but fully tabulated 
EOSs can also be implemented. As in the vacuum case, we 
employ the Komatsu-Eriguchi-Hachisu (KEH) method [76- 
81] on multiple patches [82] in order to be able to treat bi¬ 
naries consisting of different compact objects. The multiple 
coordinates systems used in this work are not necessary for 
the computation of initial data coming from equal-mass bi¬ 
naries with spins either aligned or antialigned. In this case, 
in fact, it is possible to use a coordinate system positioned at 
the center of one compact object and employ a tt symmetry 
to acquire the complete solution. Nevertheless, we here use 
three different coordinate systems (i.e., we solve all equations 
separately in three patches) as a first step toward solving for 
asymmetric binary systems, which will be presented in a fu¬ 
ture work. The gravitational equations are solved using the 
COCAL Poisson solvers (with appropriate Green’s functions), 
while for the conservation of rest mass, we follow [83] and 
employ the least-squared algorithm demonstrating the versa¬ 
tility of the methods used by our code. In this way we can 
calculate sequences of corotating, irrotational, as well as spin¬ 
ning binaries, where for the last case we use the formulation 
of [63] after small modifications to adapt it to our numerical 
methods. 

The paper is organized as follows: In Sec. II we discuss 


the equations to be solved and the assumptions made, both 
for the gravitational field in Sec. II A, as well as for the fluid 
part in Sec. IIB. Eor the latter we present the forms used in 
COCAL code for corotating, irrotational, and spinning cases. 
Section III represents the core of this work. In Sec. Ill A we 
briefly review the gravitational multipatch coordinate systems 
used in [73, 74] and discuss additional changes that are related 
to the neutron-star surface. In Sec. IIIB we describe the re¬ 
moval of dimensions from the equations, in conjunction with 
the scaling introduced in Sec. Ill A. Section III C describes the 
Green’s function used for the star patch, while in Sec. HID 
the least-squared method for solving a spinning configuration 
is introduced. Tests for our new code are presented in Sec. 
IV A for corotating binaries and in Sec. IV B for irrotational 
ones, while spinning solutions with piecewise polytropes are 
presented in Sec. IV C. A number of appendixes provides 
more technical details on several topics. More specifically. 
Appendix A reports the expressions used for the calculation 
of the mass and angular momentum of the binary. Appendix 
B illustrates a different approach to obtain a solution of the 
Tolmann-Oppenheimer-Volkoff (TOV) equations, while Ap¬ 
pendix C describes in detail the full iteration scheme and Ap¬ 
pendix D shows tests of COCAL in a very different regime 
of compactness by considering binaries of white dwarfs. Ei- 
nally. Appendix E reports the post-Newtonian expressions for 
the binding energy and orbital angular momentum of a binary 
system in quasicircular orbit, which are used as a reference. 

Hereafter, spacetime indices will be indicated with Greek 
letters, and spatial indices with Latin lowercase letters. The 
metric has signature —h -f-l-, and we use a set of geometric 
units in which G = c = Mq = 1, unless stated otherwise. 

II. QUASIEQUILIBRIUM EQUATIONS 

In this section we review the basic equations that need to be 
solved to obtain binary equilibrium configurations. Details of 
the initial-data formalism can be found in [84-87]. Here, we 
only mention the points that are relevant to the COCAL ’s new 
developments. 

A. The gravitational equations 

One of the most fruitful ideas in simulating the circular mo¬ 
tion of two bodies in general relativity was the introduction 
of helical-symmetry approximation [88, 89]. Solutions with 
such symmetry are stationary in the corotating frame and have 
a long history, starting from the electromagnetic two-body 
problem [90]. Analogous solutions in the post-Minkowski ap¬ 
proximation have been derived in [91, 92] . Helical symmetry 
was also used to obtain the first law of binary star thermody¬ 
namics [93], as well as to produce equilibrium configurations 
of binary black holes [94, 95]. 

Neglecting the loss of energy due to gravitational radiation 
and assuming closed orbits for the binary system, results in 
the existence of a helical Killing vector 

kt^ := , 


( 1 ) 
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such that 


^kgaP = 0 , ( 2 ) 

where is the Lie derivative along k and is the genera¬ 
tor of rotational symmetry. In a Cartesian coordinate system, 
but without loss of generality, we can assume the generator of 
rotational symmetry to have components 


= (-y, a:, 0). (3) 

Writing the spacetime metric in 3 h- 1 form as [37, 87, 96, 97] 

ds^ = —a^dt^ + jij{dx’' + j3^dt){dx^ + /3^dt ), (4) 

where a, /3®, 7 ^^ are, respectively, the lapse, the shift vector, 
and the three-metric on Sj, the generator of time translations 
in the rotating frame can be expressed as 

:= an^ + . (5) 

Here, uj^ := -f is the corotating shift, and the unit 
normal to Et, := Since ^klij = 0 = ^kKij, 

the evolution equation of 7 ^^ specifies the extrinsic curvature 
in terms of the shift and the lapse 

Kij — ^ DjUJi^ , (6) 

where D the derivative operator associated with 7 ^. The 

traceless extrinsic curvature 

Aij := Kij — —K ^'jij = Kij — j (7) 

can be written in terms of the longitudinal operator L 

If 2 f, 

Aij = ^ \DiUJj DjiOi —'jijDkCj 

:= — (Lw)ij . ( 8 ) 

Assuming a maximal and conformally flat slice [38, 39] 

lij='ip^dij, (9) 

the traceless extrinsic curvature Eq. ( 8 ) is written in terms of 
the shift 

2a \ 3 

= V(L/3)*^. ( 10 ) 

2a 

The tilde symbol on the longitudinal operator L denotes the 
fact that it is related to the conformally flat geometry. In de¬ 
riving Eq. (10) use has been made of the fact that dif)^ = 
0 = df^i(pj') [cf., Eq. (3)]. In the conformally flat geometry, 
the contravariant components of the shift remain the same as 
in the original spatial geometry (i.e., = /?*), while this is 

not true for the covariant components. 

With the help of Eq. (10), the constraint equations and the 
spatial trace of the time derivative of the extrinsic curvature 


result in five elliptic equations for the conformal factor ip, the 
shift /3®, and the lapse function a 

( 11 ) 

7iIA ~ ~ 

V'M) = ^{hpr^{hpy^,aSjb + 27raiP^{E + 2S) 


■.= s^ + sl, 

v 7 ‘ = 

:=S| + S'. 


( 12 ) 

(L/3)*-^ -|- 167ra'0"^j* 

(13) 


We have denoted with Sf the sources of the Poisson-type 
equations that come from the energy-momentum tensor, while 
with S'® are the sources that come from the nonlinear part of 
the Einstein tensor Gajs- The matter sources in Eqs. (11)- 
(13),S^,S/, S^, are related to the corresponding projections 
of the energy-momentum tensor 

E := , (14) 

S := 7a/3T“'^ , (15) 

f := , (16) 


where E is the energy density as measured by a “normal” ob¬ 
server, that is, an observer with four-velocity n. Note also that 
since we use G = c = Mq = 1 , all quantities in equations 
(11)-(13) are dimensionless. This is contrary to some previ¬ 
ous works (e.g., [46]), where only G = c = 1 was assumed, 
and a procedure to remove units was applied through the use 
of the adiabatic constant K. Here the normalization scheme 
used is explained in detail in Sec. IIIB. 

The above set of equations must be supplied with condi¬ 
tions on the boundary of our computational region. Since we 
will consider only binary stars, our boundary for the gravi¬ 
tational equations will be only at spatial infinity, where we 
impose asymptotic flatness, i.e.. 


lim '0 = 1, lim a = 1, lim /3* = 0 . (17) 

r—¥oo r—^00 r—^co 


We recall that a helically symmetric spacetime cannot be 
asymptotically flat, because a helically symmetric binary pro¬ 
duces an infinite amount of radiation. Therefore conditions 
(17) seem to contradict assumption (1). In reality, the heli¬ 
cal symmetry is only an approximation that is valid either for 
long times when the binary is widely separated, or for only a 
short time when the binary is tight. In practice, the emission 
of gravitational radiation reaction will lead to an inspiral, thus 
breaking the symmetry. 


B. The fluid equations 

Let be the four-velocity of the fluid. We consider a per¬ 
fect fluid with energy-momentum tensor [37, 98] 

Tai3 = (e + p)uaUp + pgap = phu^up + pgap , (18) 
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where p, e, h, and p are, respectively, the rest-mass density, the 
total energy density, the specific enthalpy, and the pressure as 
measured in the rest frame of the fluid. The specific internal 
energy e is related to the enthalpy through' 


h:== l + e+-. (19) 

P P 

The first law of thermodynamics, de = pTds + hdp, where 
s is the specific entropy, written in terms of the specific en¬ 
thalpy h, reads dh = Tds + dp/p. For isentropic config¬ 
urations, like the quasiequilibrium solutions we are seeking, 
given an EOS which relates, for example, the pressure p with 
the rest-mass energy density p, we can see that the extra vari¬ 
ables that enter our problem are p (or p) and the four-velocity 
u°‘. For these new variables, extra equations need to be used 
exploiting conservation laws. In particular, from the conser¬ 
vation of the energy-momentum tensor 

0 = = p[u°‘V a.[hu^) + 

+hu^Va{pu^) - pTV^s , ( 20 ) 

assuming isentropic configurations and local conservation of 

rest mass^ 

V„(pu“)=0, (21) 

we arrive at the relativistic Euler equation 

u°‘Va{hup) + Vph = 0. ( 22 ) 

Although we use four-dimensional indices, this is a fully spa¬ 
tial equation, since the projection along the fluid flow is triv¬ 
ially satisfied. Equations (21) and (22) provide us with four 
more equations for the fluid variables. If one of them is used 
for the determination of p, we are left with three equations that 
must determine the four-velocity m“, which has only three in¬ 
dependent components. 

Expressing the four-velocity as = m*(1, v'^) and in anal¬ 
ogy with a Newtonian decomposition, we can split the spatial 
component u® into two parts: one that follows the orbital path 
(fp, and one that represents the velocity in the corotating frame 
E*. Using the helical Killing vector, Eq. (1), we can write 

= (23) 

where the spatial part, E“ = ( 0 , E*) can be considered to be 
the “nonrotating” part of the fluid flow (i.e., the velocity in the 
corotating frame). The conservation of rest mass (21), and the 


^ Note that many authors use e and e to indicate the energy density and the 
specific internal energy, respectively. 

^ More precisely, the Euler equation is the projection orthogonal to the fluid 
flow of the conservation of the energy-momentum tensor and leads to three 
distinct spatial equations. On the other hand, the projection along the flow 
of = 0 yields a single equation expressing energy conservation 

[37], 


Spatial projection of the Euler equation (22), written in 3+1 
form translate to 

■^k{pu^) + -D,{apu^v^) = 0, (24) 

a 

+ A + hujV^^ 

+V^{Dj{hui) - D,{huj)) = 0. (25) 

The last term of Eq. (25) involves the relativistic vorticity 
tensor [37] 

Wa/3 := Vaihup) - Wp{hUa) (26) 

and is zero for an irrotational flow. It is not difficult to show 
that in the presence of a generic Killing vector field (e.g., the 
helical Killing field) k, the following identity holds [37] 

^u{hu-k) = 0. (27) 

In the case of a rigid corotation of the binary system, u = 
k, so that the Lie derivative along the fluid four-velocity u 
in Eq. (27) can be replaced by the Lie derivative along the 
helical Killing vector k. This yields J^'kh = 0 and expresses 
that in the corotating frame the fluid properties do not change. 
When the fluid four-velocity does not coincide with the helical 
Killing field, but the two vector fields are not too different, 
i.e., when u, ~ fc, expression (27) can still be true and indeed 
for the flows we will consider hereafter we will assume the 
following assumption 

= 0 = ^k{pu^) ■ (28) 

While Eq. (28) is an assumption, its correctness can only be 
assessed a posteriori and could indeed not represent a valid 
approximation if the stars are spinning very rapidly, a case we 
will not investigate here. 

In the following we will specialize Eqs. (24) and (25) under 
the assumptions (28) for corotating, irrotational and slowly 
rotating flows. Before closing we introduce a quantity that 
will be used often in subsequent sections, namely, the spatial 
projection of the specific enthalpy current 

Ui := ■ (29) 

1. Corotating binaries 

The corotating case, also called of rigid rotation [99, 100], 
is the simplest case, since the spatial fluid velocity E® van¬ 
ishes, = u*k°‘, and thus the fluid is at rest in a corotating 
frame. This means that apart from the gravitational variables 
we have only two extra fluid variables, for exam¬ 
ple p and u* once an EOS is fixed. The conservation of rest 
mass (24) is trivially satisfied, while the Euler equation (25) 
becomes a single integral equation that, together with the nor¬ 
malization condition u°‘ua = — 1, will determine all our fluid 
variables. 

In particular, the specific enthalpy current becomes 

id = hu*uj'^ 


(30) 
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and from the four-velocity normalization condition we have 


u^ = 


'Jo? — 




(31) 


Equation (25), on the other hand, has the first integral 


4 = C, (32) 

where C is a constant to be determined. Equations (31) and 
(32), together with the gravitational potentials, completely de¬ 
termine the solution for this case. We note that in all equa¬ 
tions to be solved (the gravitational ones included) two con¬ 
stants are involved. One is C, the constant that comes from 
the Euler integral, and one is 0, the orbital angular velocity. 
Thus, in order to be able to achieve a solution for our system, a 
self-consistent scheme that involves the determination of both 
Cand O must be employed. As we will elaborate later on, this 
will be achieved in conjunction with the determination of the 
length scale i?o of our problem. 

Eor the corotating case the matter sources in Eqs. (11)-(13) 
are 


E = — q\ , (33) 

E + 2S = p[h[i{au^f - 2] -b 5q \, (34) 

f = pau^iE , (35) 

where q := p/p. As we have already mentioned, all quantities 
appearing above are dimensionless, while in previous studies, 
where geometric units were used, Eqs. (33)-(35) had units of 
length”^. 


2. Irrotational and spinning binaries 

Irrotational configurations have ujap = 0 , so that the spe¬ 
cific enthalpy current hua can be derived from a potential [37] 
i.e.. 


hUa = 'S/= Da^ + Ua^n^ , (36) 

SO that Ui = Di^ since 7 • n = 0. To allow for spinning 
configurations we need to extend expression (36) and we do 
this following Ref. [67] and introducing a four-vector s (not 
to be confused with the specific entropy s) 

hUa = Vq,*!) -b Sa , (37) 

so that M* is decomposed as 

+ s*, (38) 

where the part corresponds to the “irrotational part” of 
the flow and the s* part to the “spinning part” of the flow. In 
what follows we will present expressions for spinning binaries 
(s® J 0 ) and one can recover the irrotational ones by setting 
the spinning component s* equal to zero. 

Using the decomposition (38) and the assumption (28), the 
Euler equation (25) can be rewritten as 

!L + VW,<^]=0, (39) 

J 


and hereafter we will assume 

= 0 , (40) 

which is likely to be a very good approximation in the case 
of slowly and uniformly rotating stars, for which Si is intrin¬ 
sically small. ^ Hence, the Euler equation for generic binaries 
(25) 


7 “ [^k{huoJ + JSfv(sa)]+A = 0, (41) 

under the assumptions (28)i and (40), yields the reduced Euler 
integral 


+ = C, (42) 

where again C is a constant to be determined. A few remarks 
should be made at this point. Eirst, it is not difficult to obtain 
the following identity 

+^v{Sa)\ = 

Tz a^) “b (^ck) ~b ; ('^3) 

SO that our assumptions (28) i and (40) 

{hua) = 0 = 7 “.if^ (so), (44) 

are equivalent to setting the left-hand side of Eq. (43) to zero. 
In turn, this implies that also the right-hand side of (43) is 
zero, which is true if, for instance, each of the three terms is 
zero, i.e., if 

7“,5ffc(Va«') = 0 = 7“-^V<[./(h«*)(Sa) = lf^s/{hu*){Sa) ■ 

(45) 

The three conditions in (45) coincide with the assumptions 
made in [63]. Stated differently, because the conditions (44) 
are compatible with the conditions (45), it does not come as a 
surprise that we obtain the same Euler integral (42) as in [63] 
despite making apparently different assumptions [cf., (44) vs 
(45)]. Second, using the decomposition (37), it follows that 

-if^k{huoJ = 7 “[^fc(V„$) + ^k{soJ ], (46) 

and hence the question about EEk{hua) = 0 depends on both 
-2fc(Vo,4)) and JEk{sa) being zero^. The second term is es¬ 
sentially an input to our problem, while the first one comes 
from the conservation of rest mass (57), which depends on the 
spin input Sa- Einally, although the Euler integral has the same 
form for both irrotational and spinning binaries, it produces a 
different equation since the three-velocity U® is different in 
these two cases. More specifically, it is 

= hu\uj^ + V^), (47) 


^ In practice we will consider stars with spin period down to 0.6 ms, but 
this is still ’’slowly” spinning when compai'ed to the minimum period. 

^ Note that even when the spins are aligned with the orbital angular momen¬ 
tum ^^( 5 ^) ^ 0. 


+ Di 
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so that 


= 




hu* 


— OJ 


(48) 


Neumann type, that is, in terms of derivatives of the rest-mass 
density and of $ 

= 0. (58) 


In this case, the fluid variables are p (or equivalently p or 
h), u*, and the fluid potential fl). The equations that will deter¬ 
mine them are the normalization condition = — 1, the 

Euler integral (42) [with the use of Eq. (48)], and the conser¬ 
vation of rest mass (24). 

In particular, from the norm of ft* we get 

h = ^Ja'^{hu^Y - (A4> + s*)(E>*4> + s*), (49) 


therefore, the Euler integral (42) takes the following form 
quadratic in hu* 

a^{hu*f - A(hu*) - -p s*) = 0, (50) 

where \ = C + uj'^Di^. Thus 


hu* = 


A + + 4a2si(i9®<l> -I- s*) 


(51) 


where we take the positive root since the negative one is in¬ 
correct at least in the limit of s® = 0, when it yields hit* = 0. 

Having computed A, we first calculate hu* from Eq. (51), 
and then h from Eq. (49). Eor purely irrotational binaries 
hu* = X/a^ and h = ja^ — 

Although we will not make immediate use of u* and h sep¬ 
arately, we report below their form for completeness 

h= -{D,^ + Si){D^^ + s^), (52) 

„,t _ + (A4> + -f s*) 

ha 

where 


:= 


A^ -f 2a^W -f AVA2 -f 


IE := -f . 


(54) 

(55) 


The potential <1> will be computed from the conservation of 
rest mass (24), which under Eq. (48), and after expressing the 
spin velocity as a power law [63] 


s* = , A G K (56) 


will produce an extra elliptic equation 


W 


- 4 ’ 


A+A 


dii 


*9An 


a pip' 


.&+A 


= S^. 


(57) 


The boundary for the fluid is represented by the surface of the 
star; hence the boundary condition for Eq. (57) will be of von 


A possible and convenient choice for the parameter A that 
will be used in Sec. IVC is A = —6, as it removes the last 
term in Eq. (57)^. Any other value will not change the charac¬ 
ter of the equation or the boundary condition, although it will 
change the detailed properties of the flow velocity and there¬ 
fore of the binary. We will comment on this point in Sec. IV C, 
where we will also illustrate the results for A = 0. 

Eor the spinning case, the matter sources in Eqs. (11)-(13) 
are 


E = p 
E + 2S = p 


ra 2 


{huy-q 


3 — (huy — 2h + 5q 
h 


(59) 

(60) 


f = pau^vX = p'^ihvP) \4 ^19*$ + yy , (61) 


where we have used Eq. (51) to simplify the calculations. 


C. Equation of state 

The EOS used in this work is represented by a sequence 
of poly tropes called a piecewise polytrope. This is proven to 
be a good approximation for a great variety of models [101- 
104]. If N is the number of such poly tropes, in each piece the 
pressure and the rest-mass density are 

p = Ky\ t = l,2,...,N. (62) 

The order of the polytropes is i = 1 for the crust, and i = N 
for the core, and Eq. (62) holds for 

Pi-1 < P < Pi- (63) 

As we have discussed in Sec. IIB, the first law of thermody¬ 
namics for isentropic configurations gives dh = dp/p, which 
can be expressed in terms of q to yield 

dh = -^^dq, (64) 

or equivalently 

h-hi= p ^ (g- gi)) (65) 

where hi,qi correspond to values at the right end point (the 
one closest to the core) of the ?-th interval. In terms of q, we 


^ More precisely, A = —6 makes the spin have zero divergence in the three- 
geometry {Dis'^ = 0) if we choose it to have zero divergence in the con¬ 
formal three-geometry {Dis'' = 0). 
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can express the rest of thermodynamic variables as 

p = , ( 66 ) 

e = ph — p. ( 68 ) 

Enforcing the continuity of the pressure at the A^ — 1 interfaces 
of each interval constraints all adiabatic constants Ki but one 

K,p^^ = . (69) 

As a result, the free parameters are; one adiabatic constant, 
— 1 rest-mass densities, and N adiabatic indices, a total of 
2N parameters. 


III. NUMERICAL METHOD 

The COCAL code for binary black holes has been described 
in detail in Refs. [73-75]. Here, we will review the most 
salient features of the grids used for the solution of the field 
equations and discuss the differences that arise from the treat¬ 
ment of binary stars. 

When treating binary systems, COCAL employs two kinds 
of coordinate systems. The first kind is compact object coor¬ 
dinate patch (COCP) and has exactly two members, one cen¬ 
tered on each star. The second kind is asymptotic region coor¬ 
dinate patch (ARCP) and can have in principle any number of 
members in an onion type of structure. In our computations, 
the ARCP patch has only one member, which is centered on 
the center of mass of the system. All coordinate systems use 
spherical coordinates G [ra,rb] x [0,7r] x [0, 27r],but 

the components of field variables (like the shift) are written 
with their Cartesian components (/?^, /3^, /3^). The values of 
Ta, Tfc depend on the compact object (black hole or neutron 
star) and the coordinate patch (COCP or ARCP). For the bi¬ 
nary systems treated here, of the COCP patch will always 
be zero, while for the ARCP patch C>(10M), M being 
the mass of the star. As for the values are kept the same as 
in the binary black hole computation, i.e., 0{\Q^M) for the 
COCP patch, and C>(10®M) for the ARCP patch. 

The orientation of the coordinate patches is as follows; the 
ARCP patch has the familiar (x, y, z) orientation, the first 
COCP patch, which is centered on the negative x axis of the 
ARCP patch, has the same orientation as the ARCP patch, 
while the second COCP patch, which is centered on the posi¬ 
tive X axis of the ARCP patch, has negative (x, y) orientation, 
and positive z orientation with respect to the ARCP patch and 
to the first COCP patch. In other words the coordinate system 
of the second COCP patch is obtained from the first COCP 
patch by a rotation through an angle of tt. 

The geometry of the ARCP patch (or any number of them) 
is that of a solid spherical shell with inner radius and outer 
radius rf,. On the other hand, the COCP patch geometry is 
that of a sphere of radius rt, with another sphere of radius re 
at distance dg from the center, being removed from its inte¬ 
rior. This second sphere whose boundary we call the excised 
sphere Se, is centered on the x axis around the other compact 


Ta'- Radial coordinate where the radial grids start. For 
the COCP patch it is Va = 0. 
ri,: Radial coordinate where the radial grids end. 

Tc'- Center of mass point. Excised sphere is located 
at 2rc in the COCP patch. 

re'- Radius of the excised sphere. Only in the COCP patch. 
rs'. Radius of the sphere bounding the star’s surface. 

It is Ts < 1. Only in COCP. 

Nr'. Number of intervals Ar; in r £ [ra,rb\- 

Nr'. Number of intervals Ar^ in r £ [0,1]. Only 
in the COCP patch. 

Nr'. Number of intervals Ar; in r £ [0, rs] in the COCP patch 
or r £ [ra , r^ -f 1] in the ARCP patch. 

A™: Number of intervals Ar; in r £ rc]. 

Ng: Number of intervals A6j in 0 £ [0, tt]. 

N^: Number of intervals A(j)k in (j) £ [0, 27r]. 
d: Coordinate distance between the center of Sa (r = 0) 
and the center of mass. 

da: Coordinate distance between the center of Sa (r = 0) 
and the center of Se- 
L'. Order of included multipoles. 

TABLE I. Summary of grid the parameters used for the binary sys¬ 
tems computed here. 


object. For the first COCP patch, its excised sphere Se is cen¬ 
tered on the position of the second star, while for the second 
COCP patch, its excised sphere Se is centered on the position 
of the first star. The size of every sphere Se is slightly larger 
than the star resolved with an opening half-angle of ~ 7 r /3 
as seen from the origin of the COCP patch. This is done to 
resolve accurately the region around the other star and reduce 
the number of multipoles used in the Legendre expansion. In 
this work we typically use 12 multipoles in our computations. 
Table I summarizes the properties of the various coordinate 
patches used and which are illustrated schematically in Fig. 
1 . 

On the coordinate grids described above we solve the 
Poisson-type nonlinear equations (11), (12), and (13). These 
equations are solved using the representation formula with a 
suitable choice of a Green’s function (this is also known as 
the KEH method [81]). The Green’s function is expanded in 
terms of spherical harmonics and the integrals (both volume 
and surface) are performed on the spherical coordinate grids 
described above and in more detail in Sec. Ill A for the case 
of binary stars. 

The numerical approximations introduced in COCAL are of 
different types; First, from the truncation of the series of the 
Legendre expansion; second, from the solution of the equa¬ 
tions on discretized grids via finite-difference methods [37]. 
Typically, we use second-order centered stencils for the nu¬ 
merical differentiations and integrations [73]. An exception 
is the use of a third-order finite-difference stencil for the ra¬ 
dial derivatives and the use of a fourth-order integration in the 
polar coordinate [74]. We also typically use second-order in- 
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terpolations of scalar functions from grid points to midpoints. 
Furthermore, when a function needs to be interpolated from 
one coordinate system to another, we use a fourth-order La¬ 
grange interpolation. This usually happens when we compute 
the surface integral at the excised sphere Se for in that case 
we need the potential and its derivative on Se as seen from the 
center of Se- 


A. The numerical grids 

In Refs. [73, 74] we described in detail optimal numeri¬ 
cal grids that we constructed in order to lower the error of the 
potentials for both close and large separations, for any kind 
of mass ratio. In particular, when we computed sequences, 
instead of keeping the radii of the black holes the same and 
increasing their separation, we kept the separation fixed and 
decreased their radii. By choosing the interval separation near 
the black holes according to their excised radius, we were able 
to obtain sequences comparable to the ones produced by spec¬ 
tral methods. 

We adopt here the same philosophy for the computation of 
binary stars. Contrary to previous studies [46], in order to 
compute sequences of binary stars we let the maximum radius 
of the star be variable and we denote by Vg the infimum of 
the radii of all spheres bounding the star that are centered on 
the origin of the COCP patch. By continuously diminishing 
Tg while keeping the distance between the stars constant, we 
can compute sequences of stars with a continuously increasing 
separation. In this way we can control the region around the 
excised sphere as described in [73, 74], while maintaining the 
accuracy in the area covered by the neutron star. 

The COCAL radial grid for binary stars can be seen in Fig. 
1, where all the radial distances are the normalized quanti¬ 
ties discussed in Sec. IIIB. In that sense, they should be de¬ 
noted with a hat, for example f, which is omitted here for 
simplicity. When comparing with Fig. 2 of [73], we can ob¬ 
serve that there is also an important difference in the notation. 
This regards the quantity N^, which previously was used to 
denote the number of intervals in 1], while here is used 
to denote the number of intervals in [ra,rs] = [0,r,j], with 
Tg < 1. The number of intervals in [0,1] is denoted by a 
new variable called (this plays the role of the old for 
grid comparisons). This change was necessary since in all 
previous studies [46], the surface of the star was bounded by 
the fixed r = 1 sphere, and therefore the fluid extended until 
that point. To compute stars at larger separation while satis¬ 
fying this constraint, we would have to increase the position 
of the excised sphere Se and therefore expand all grid quanti¬ 
ties analogously. To avoid such a complication, we introduce 
a variable Vg that effectively mimics the change in separation. 
By varying the end point of the fluid (point Vg) we achieve 
the same result as varying the distance between the stars, but 
we maintain the good convergence properties that were es¬ 
tablished in [74] while maintaining our fluid code essentially 
unchanged. 

As we can see from Fig. 1, there are five regions in the 
COCP patch of a star that are denoted by S, I, II, III, and 


IV. The star is resolved by a constant grid spacing Ar = 
fg/Nl, as region II, which has spacing Ar 2 = I/N^. Setting 
Ari := n — Ti-i, the grid intervals in each of them are 

Ari = Ar, for i = 1, • • • , A,'! - 1, (70) 

Ari+i = hi An , for i = N^, ■■■ ,N^ - 1, (71) 

Ar, = Ars, for i = a;,--- ,A“, (72) 

Ar,+i = hsAn , for i = A“, • • • , A“ + A^ - 1, (73) 

Ar,+i = hiAn , for i = A“ + Ai, • • • , A, - 1, (74) 


which correspond to regions S, I, II, III, and IV, respec¬ 
tively. The ratios hi{> 1) {i = 1,3,4) are, respectively, de¬ 
termined from the relations 


1 — r, = Ar 


2rr = Ar 


n - ire = Ar 


'^-1) 
hi-I 

h,{hf - 1 ) 

hs-l ’ 

, ,,Nr-N^-Nl 

h/i{h^ — 1) 

Ii4 — 1 


(75) 

(76) 

(77) 


For the ARCP coordinate system, there are in general two 
regions, one with constant grid spacing and one with increas¬ 
ing spacing. The grid intervals in these regions are defined 
by 


Ari= Ari, for i = 1, • • • , A“ , (78) 

Ari+i=fcAri, for i = A“, • • • , A,. - 1, (79) 

where Ari = 1 /A);, and the ratio k is determined from 

, - 1) 

n - A =: Ar—^-. (80) 

As regards the angular resolution, we keep the same grid 
interval in the 6 and (p directions and therefore 

A9, = e, - e,.i = A0 = ^, (81) 

27r 

= (t)k- (t>k-i = A(j) = — . (82) 

One of the additional complications of having to deal with 
the fluid of a star, instead of a vacuum spacetime, is the need 
to accurately And its surface. This surface may contract or 
expand during the calculation, creating significant problems 
in close binary configurations. One very effective solution 
to these issues [105] is the use of surface-fitted coordinates 
(SFC) that exist only inside each fluid and are normalized by 
the radius of the star. We denote this extra spherical coordinate 
system as {rf. Of, fif), where 

and where the surface of the star is denoted by R{0,(j)). 

By construction, the domain of these fluid coordinates is 
[0,1] X [0, tt] X [0, 27r], and R{9f,(j)f) is a function that will be 
determined at the end of the self-consistent iterative method 
(see Appendix C). The advantage of SFC in the computation 
of derivatives on the star’s surface, as well as the implementa¬ 
tion of the boundary condition Eq. (58), will be discussed in 
Sec. HID. 
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S'. Equidistant, intervals. 

/: Non-eqnidistant. Intervals in 5 + /: intervals, 

11 Equidistant. Intervals in 5 + / + //: N™. 

Ill: Non-eqnidistant, intervals. 

IV: Non-eqnidistant, Nr — {Nr + Ai™) intervals. 


FIG. 1. Structure of the radial grid for the COCP coordinate system for binary stars, can take any value in (0,1], but typical values used are 
in the range [0.5, 1]. Decreasing Ts amounts to a larger effective distance between the two stars. 


B. Dimensionless and normalized variables 

Having removed the dimensions from our equations by us¬ 
ing units in which G = c = Mq = 1, we perform a normal¬ 
ization of all variables in order to introduce a scale in our prob¬ 
lem that is intimately related to the variable Vg, introduced in 
Sec. Ill A. 

We normalize variables by demanding that the intersection 
of the star’s surface with the positive x axis be at r = Vg. If 
Rq is the scaling parameter, we impose® 


fi(V2,0) 

Rf) 


^(7r/2,7r) 

Ro 


(84) 


so that VgRo is the real semimajor radius of the star. Hereafter, 
we will denote normalized variables with a hat and thus define 


X := 




(85) 


0 


from which it follows that the normalized version of Eq. (11) 

is^ 


V^i; = S^.+RlSi, 


( 86 ) 


® Note that rg, but also rj in Eq. (83), are ratios of two radial coordinates 
and thus dimensionless for any choice of units; in this respect, they do not 
need to be indicated with a hat. 

’’ Similar noiuialized equations hold for Eqs. (12) and (13). 


where V is the Laplacian operator associated with the vari¬ 
ables x^, and similarly has all derivatives with respect to 
the normalized variables. Note also that 


W* = /3* -F = /3* + , 


where Cl := flRn. If we now define 


$ := 


$ 

Ro 


(87) 


( 88 ) 


we observe that all scaling factors in Eq. (57) drop out. Be¬ 
cause this is also true for the boundary condition, i.e., Eq. 
(58), these equations are each coded in the same form but with 
normalized quantities replacing the original ones. 

Before proceeding further with our normalization scheme, 
let us comment that the surface-fitted coordinates of Eq. (83) 
are already normalized coordinates and their radial range 
is [0,1] irrespective of the fluid scaling profile Xg. Since 
R{S,(j)) < i?(7r/2,0) = RoXg, we have 0 < r < R[9,(j)) < 
RoXg, so that the range for f is 


r < R{0, cf)) <rg, 


Tf = 


R{0,C) 


< 1 ■ (89) 


Changing the scale by modifying Rq will affect the conformal 
factor and the lapse function since they scale as [46, 106] 


ip = 


a = a^° 


(90) 
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As mentioned earlier, the system of partial differential 
equations that we have to solve i.e., the normalized versions of 
Eqs. (11)-(13) and Eq. (57), involve three constants: Rq, Ct, 
and C. To find them we will use the Euler integral evaluated 
at three arbitrary points to construct a nonlinear 3x3 system 
that will be solved with a typical Newton-Raphson method. 
This procedure will be repeated every time we solve for any 
of the unknown variables ■0, /?*, a, $, and q, since any change 
of them will affect the Euler integral and thus the three con¬ 
stants. 

The arbitrary points we choose to evaluate the Euler integral 
are along the x axis and in spherical coordinates are defined 
as 


ri 

= Ts , 

9i 

= 77 / 2 , 

</>i 

7-2 

= 0, 

92 

= 0, 

02 

rs 

= rs , 

93 

= 77 / 2 , 

03 


(91) 

(92) 

(93) 


In the corotating case, Eqs. (32) and (31) will become 


Fc{Cl, Rq, C) = — InC + i?o Ind + In/i 


+ iln 


02 


2Ri 




1- ^ (/3" + ^20") 


= 0, (94) 


while for the spinning case, Eqs. (42), (49), and (51) yield 


Fis{Cl, Rq, C) = — In a + i?Q Ind -f In /i 

1 


— In 


1 


+ 


AiRo) 

A2 


■In 


1 + 


B{Ro) 


= 0(95) 


with 

A{Rq) := 5 * 4 $ + 

B(i?o) := + 2(0^)^' 04$ , 

X ■= c + {py + vt^y)dy^. 


Evaluating either Eq. (94) or Eq. (95) at the three points 
given by (91)-(93) will produce a system of three equations 
in the three unknowns Rq, and C; the solution of the 
system will determine these constants. These unknowns are 
computed separately for the two stars, although they yield the 
same solution for the case of equal-mass binaries considered 
here. We also note that the star’s surface remains fixed and the 
values of the specific enthalpy h at the three points (91)-(93) 
is unchanged, with /i = 1 at (91) and (93). 


C. Elliptic solver for the gravitational part 

As discussed in detail in Refs. [73-75], Eqs. (11)-(13) are 
solved using the representation theorem of partial differential 
equations in a self-consistent way. Starting from 


( 96 ) 


where S' is a nonlinear function of /, and using the Green’s 
function without boundary G(x, x') = l/|x — x'| that satisfies 

V^G{x,x') =—'iTT5{x — x'), (97) 

a solution for / is obtained as 

/(x) = — — [ G{x,x')S{x')(f’x 

47r Jy 

+ ± f [G(x, x')V'“/( 0 ) - /(a:')V'“G(x, x')] dS', . 

47^ Jdv 

(98) 


where V is the domain of integration, x,x' G V Q Eq, the 
initial spacelike hypersurface. The volume V and its bound¬ 
ary dV depend on the coordinate system we are solving for, 
i.e., COCP or ARCP as described in Sec. Ill A. This is the 
principle of the KEH method [81] and will be suitably mod¬ 
ified in order to account for the specific boundary conditions 
that exist in the new COCAL coordinate systems. Eor exam¬ 
ple, the conformal factor will be expressed as 


0(x) = x{x) + 0iNT(a:), (99) 


where 


0INT(x) = -^j^. 


1 f Sl{x') + S^^{x') , 


1 

-b-— 


477 JdV I \x-x 




d'^x' 

1 


X — X' 


dS',, (100) 


and 


X{x) = ^ [ [G®°(x,x')V'“(0bc - 0int(0) 

477 JSaUSt 

-(0BC - 0int)(0)V'“G^'^(x,x')] dS'^ . (101) 

Note that G®^ is the Green’s function associated with the 
boundary conditions applied on the corresponding field 0 bc 
at the boundaries Sa and Sb- Eormulas (99)-(101) will be ap¬ 
plied separately on every coordinate patch. If, for example, 
we have one ARCP patch (as it happens in our computations) 
it means that the equations above will be applied three times: 
two for the COCP patches and one for the ARCP patch. Of 
course, the domains of integration vary according to the dif¬ 
ferent patches considered. More specifically, if we denote by 
B{R) a sphere of radius R in each of the COCP patch, then the 
integration domain of Eq. (100) will heV = B{rb) — B(re) 
and dV = S'e U S't, = dB{re) U dB{rb), while that of Eq. 
(101) will be SaU Sb = Sb, since = 0 for star configu¬ 
rations in the COCP patch. Similarly, in the ARCP patch the 
integration domain of Eq. (100) will be F = B{rb) — Blja), 
and that of Eq. (101) dV = Sa^ Sb- 

We recall that in Ref. [73] we have introduced a number of 
Green’s functions G^^{x,x') suitable for various boundary 
conditions. Here we add one more Green’s function used in 
the COCP patch 


GSD(x,x') :=^gfD(r,r')^e, 


^=0 


m—0 


(£ — m)! 


xP(^{cos9) P^'^{cos9’)cos[m{(j) — 0 ')] ,( 102 ) 
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where 


9r{r,r') 







(103) 


and eo = 1, em = 2 for m > 1, while are the asso¬ 
ciated Legendre polynomials, and r> := sup{r, r'}, := 

inf{r, r'}. 

In the ARCP patch we use a Green’s function x'), 

whose radial part satisfies the Dirichlet-Dirichlet boundary 
conditions on Sa and Sb 






t 




.(104) 


D. Elliptic solver for the fluid part 

Next, we describe the method used to solve Eq. (57) and 
which is therefore valid only for the spinning binaries. The 
boundary condition for $, Eq. (58), is of von Neumann type 
and therefore we could apply the Poisson solver of Sec. IIIC 
to obtain a solution. Instead, and as a demonstration of the 
versatility of the methods employed by COCAL , we will adapt 
the procedure discussed in Ref. [46] and solve this boundary- 
value problem as an application of the least-squares algorithm. 

Eirst, we assume that the solution of Eq. (57) can be written 
in the form 

= -7- f + C{x) = $y(a;) -f ({x ), 

Jv \x — X \ 

(105) 

with obeying the Laplace equation 

V^C(a:)=0. (106) 

Using the decomposition of Eq. ( 105), the boundary condition 
(58) is written as 

— m’'di^v — = m^di( , (107) 

where we used the normal to the surface unit vector = 
{nsY instead of the gradient of the rest-mass density. The 
equation above is evaluated on the surface of the star, R{9, </>), 
and the velocity potential satisfies the following symmetries 

^{r,Tr - e,(j)) = 4>(r, 0, ())), (108) 

$(r, 0, 27r - (109) 

which in turn imply that the homogeneous solution, which is 
regular at the stellar center, can be expanded as 

L I 

C,{r,e,4>) = X! X! atrnr^[l+{-lY^'^]Yl^{cose)sin{m4 >), 

m—1 

( 110 ) 

where are the spherical harmonics and aem coefficients 
to be determined. When the spins of the stars are in arbitrary 


directions, the symmetries (108) and (109) no longer apply 
and expression (110) will contain also cosine terms. 

On the other hand, if R{9,(j)) is the surface of the star, the 
spatial vector connecting any point on it with the center of 
coordinates is given by 


x{9,(j)) = R{9,(l)) (sin 9 cos (p, sin 9 sin cj), cos 9). 
thus the unit normal vector will be 




or equivalently 
1 


y/h 


1 OR. 


ris = —:= \ r - ——9 - 


dx dx 
89 dp 


dR 


Rd9 RsRv9 dp 


4> 


where r,0,p are the spherical unit vectors and 


h{9^ p) :— 1 -|- 


R~^ 


1 dR 


( 111 ) 


( 112 ) 


(113) 


(114) 


Rsm9 dp ^ 

Using Eqs. (110) and (113), the boundary condition (107) 
is written as 

L I 

EE ae^Fe^{9,p)=H{9,p), (115) 

l—l m—1 

where 

H{9,p) := p^hu^uj^rrii — m^di^v — , (116) 


and 


Fim{9,P) := [l + (-l)^+n X 

m^-^Yrsin{mP) - sin(m</>)- 

no 

— ^^YPmcos{mp) . (117) 

dp sin 9 

To solve for the coefficients aim, we consider the functional 


^-E EE aimFr{9uPj)-H{9,,pM =0, 

9i,,<pj 1—1 m—1 

(118) 

of the discretized version of the boundary condition (115), and 
demand that for fixed indices p and q 

dS 


dar] 


= 0 . 


(119) 


The minimizing condition Eq. (119) then yields 


^ ^ aim 

l,m 


X] (di , Pj )Fpg {9i , Pj ) 

'i’J 


( 120 ) 


= Y,H{9,,Pj)Fpg{9,,p,), (121) 


which is a linear system in terms of the aim coefficients. 
Eor L even, the dimensions of the system is M x M with 
M = L{L + l)/2. After determining the coefficients aim, 
the solution for the velocity potential $ is obtained from Eqs. 
(105) and (110). 
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C=0.12 



C=0.18 



C=0.12 



C=0.18 



C=0.12 



C=0.18 



FIG. 2. Quasiequilibrium sequences for corotating binary neutron stars with EOS consisting of a single polytrope with F = 2. The left column 
corresponds to models with compactness C = 0.12, while the right column corresponds to C = 0.18. M and J are the total ADM mass and 
angular momentum of the system. The resolutions used are those given in Table II. A comparison is made with the results of Ref. [1], where 
similar solutions were obtained from the spectral code LORENE [112]. 


IV. NUMERICAL RESULTS 

In what follows we report tests of our new code against 
previous results obtained by other groups and then present 
some new ones. In particular, we focus on the construction 
of quasiequilibrium sequences for corotating, irrotational, and 
spinning binaries, and produce a binary white dwarf solution 


in order to explore the weak-field limit of our code. We com¬ 
pute two main global-error indicators to measure the accuracy 
of our converged solutions; the first one is given by the re¬ 
lation where and are the Komar 

and Arnowitt-Deser-Misner (ADM) mass, respectively [107- 
111]. The second one is instead related to the first law of bi¬ 
nary thermodynamics = QdJ, where Q is the orbital 
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Type Patch 

Ta 

Ts 

Tb 

Vc 

Te 

a/ 

A/ 

A™ 

Nr 

As 

A^ 

L 

Hs2d COCP 

- 1 

0.0 

var 

10^ 

1.25 

1.125 

50 

64 

80 

192 

48 

48 

12 

COCP 

- 2 

0.0 

var 

10^ 

1.25 

1.125 

50 

64 

80 

192 

48 

48 

12 

ARCP 


5.0 

- 

10® 

6.25 

- 

16 

- 

20 

192 

48 

48 

12 

Hs2b COCP 

- 1 

0.0 

var 

10^ 

1.0625 

1.03125 

60 

64 

68 

192 

48 

48 

12 

COCP 

- 2 

0.0 

var 

10^ 

1.0625 

1.03125 

60 

64 

68 

192 

48 

48 

12 

ARCP 


5.0 

- 

10® 

6.25 

- 

16 

- 

20 

192 

48 

48 

12 


TABLE II. Grid parameters used for the corotating sequences of compactness 0.12, 0.18, presented in Fig. 2. Hs2d refers to solutions at large 
separations, while Hs2b refers to binaries with small separations. For the Hs2d case the separation between the two neutron stars is kept fixed 
at ds = 2rc = 2.5 while the surface of the neutron star (maximum value Vs) varies from rg = 0.5 to rg = 0.86 creating effectively binaries 
with separations from 5.0 = 2.5/0.5 to 2.91 = 2.5/0.86. Similarly in the Hs2b case dg = 2rc = 2.125, and the effective separations are 
from 2.125/0.74 = 2.87 to 2.125/0.98 = 2.17. 


frequency and J the orbital angular momentum [93]. Explicit 
definitions and computational algorithms for these quantities 
within COCAL are presented in Appendix A. 

Sequences of constant rest mass can be thought of as snap¬ 
shots of an evolutionary process that drives the two stars close 
to each other as a result of gravitational radiation reaction. 
At every instant in time, the rest mass of each star is con¬ 
served; furthermore, if the flow is irrotational, the circulation 
of the fluid velocity on any loop is also conserved (Kelvin- 
Helmholtz theorem [37]). It is possible to characterize these 
sequences via the properties of the stars, such as the compact¬ 
ness or the ADM mass, when the binary has infinite separation 
and each star is spherical. 

By varying the separation between the two stars and solv¬ 
ing each time all the relevant equations, we obtain solutions of 
a given central rest-mass density, and then another loop of so¬ 
lutions has to be invoked in order to And the particular central 
rest-mass density that yields a star with the desired rest mass. 
Typically, this is done by a Newton-Raphson method and it 
takes a maximum of ten iterations depending on the starting 
solution. 

In this way we can monitor important quantities like the 
binding energy of the system, which is defined as 

Eh ■■= - Moo , (122) 

and represents the total energy lost in gravitational waves by 
the system, since Moo is twice the ADM mass of a single 
isolated spherical star. 


A. Corotating solutions 

As mentioned in the Introduction, corotating states [106, 
113-115], i.e., states with zero angular velocity of the star 
with respect to a corotating observer, probably are not physi¬ 
cally realistic due to the low viscosity of the neutron-star mat¬ 
ter. Such solutions represent an important step in the numer¬ 
ical solution of the binary problem, since they provide key 
insights for the numerical implementation of a stable algo¬ 
rithm. In particular the surface-fitted coordinates, as well as 
the solution of the Euler integral, Eq. (32), can be thoroughly 
checked. This allows us to perform a calibration without 


having to worry about the fluid flow (i.e., to solve the equa¬ 
tion of conservation of rest mass). Corotating evolutionary 
configurations are known to exhibit a minimum in the mass 
and angular momentum versus the normalized angular veloc¬ 
ity ^ which was taken to denote the putative in¬ 

nermost stable circular orbit, beyond which the binary was 
thought to proceed rapidly toward a merger. In practice, fully 
general-relativistic simulations of inspiraling binary neutron 
stars do not show the existence of such an instability, reveal¬ 
ing instead that the inspiral and merger is a smooth process 
[2, 3]. Nevertheless, the presence of such a minimum repre¬ 
sents a useful test of numerical codes, as does the appearance 
of a familiar spike, similar to the one encountered in binary 
black-hole solutions, when plotting the binding energy versus 
the angular momentum of the binary. 

In Eig. 2 we present sequences of binary neutron stars that 
correspond to the compactness of C := = 0.12 and 

0.18, where R are the (ADM) mass and radius of 

each star when taken at infinite separation. The ADM and Ko¬ 
mar mass, as well as other quantities used in COCAL , are de¬ 
scribed in detail in Appendix A. The adiabatic index is E = 2 
and the polytropic constant is set to be iT = 1. In the vari¬ 
ous plots a comparison is made between the results obtained 
with COCAL and those presented in Ref. [1], where the same 
initial data were computed using LORENE , a pseudospectral 
code developed by the Meudon group [1]. As we can see, the 
relative difference in the results between the two codes is of 
the order of 0.05%, even when a medium resolution is used 
for COCAL . The grid structure used in these calculations is 
the one described in Table II. 

Similarly, in Eig. 3 we report the change in the central rest- 
mass density with respect to the one at infinity, which is poo = 
0.0922 for compactness C = 0.12, while it is poo = 0.1956 
for C = 0.18. Clearly, the central rest-mass density decreases 
as the binary comes closer, making the onset of an instability 
to gravitational collapse very unlikely [116]. In addition, as a 
measure of accuracy of these corotating sequences, we plot in 
Eig. 4 the relative difference 


Am := 


A^adm - 


Mg 


(123) 


as a function of the binary separation ds/rg- Note that all radii 
are here normalized to the scaling factor Rq and are therefore 
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dg/rg energy density 


FIG. 3. Relative change in the central rest-mass density for the coro- FIG. 5. Corotating sequences versus central energy density for dif- 

tating sequences in Fig. 2, shown as a function of separation. ferent separations. The TOV curve corresponds to infinite separation. 

Here we use = 1, and therefore Rq is the radius of the star. 





FIG. 4. Measure of the virial error for the corotating 

sequences in Fig. 2, as a function of separation. 


dimensionless, so that, e.g., the physical distance between the 
two neutron stars is dgRo- As we can see, even for the medium 
resolution used in these calculations the error is below 10“^. 
All of the quantities in the expression above have been ex¬ 
tracted from the ARCP patch, as integrals at infinity. We note 
that at present COCAL does not implement a unifying mesh, 
and this prevents us from calculating the virial error as ob¬ 
tained by Friedman, Uryu, and Shibata [93], since we are us¬ 
ing overlapping coordinate systems. We plan to revisit this 
issue in the future. 

Finally, in Fig. 5 we report sequences of corotating bina¬ 
ries for increasing central rest-mass density at different sep¬ 
arations, from ds/Ro = 4, down to separation in which the 
two stars are almost touching dg/Ro = 2.125. Here we use 
= 1, and therefore the radius of the neutron star is Rq. 
Clearly, for any given central energy density a larger mass 
is supported by the binary (supramassive solutions) when we 
move to closer configurations, once again excluding the onset 
of an instability to gravitational collapse to black hole [117- 
119]. 


B. Irrotational sequences 


As anticipated in the Introduction, irrotational neutron stars 
have been considered as a reasonable first approximation to 
describe the flow in binary configurations. In such a case, 
the total angular momentum is less than the corresponding 
of a corotating binary, since in each star there is a flow in 
the counterdirection with respect to the orbital motion®. This 
has two consequences. First the inspiral of irrotational bina¬ 
ries is faster than that of corotating ones or, equivalently, the 
gravitational-wave frequency is expected to increase with a 
faster rate for irrotational systems. To first order in the spins, 
the rate of inspiral is [ 120, 121] 


xL- 



(124) 


where ix is the symmetric mass ratio and L the unit angular 
momentum vector. From (124) we can see that 


(125) 


which is expected, since spinning binaries have to radiate also 
the additional angular momentum before they merge [122]. 

Second, in the light of the results obtained for binary black 
holes, where binaries with larger spins lead to increasingly 
spinning final black holes [123, 124], the irrotational binary 
system will eventually lead to a Kerr black hole that is more 
slowly rotating than the corresponding one produced by the 
corotating binary. 


* For this reason these binaries are also called counter-rotating configura¬ 
tions. 
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FIG. 6. Sequences of irrotational binary neutron stars, with an EOS consisting of a single polytrope with F = 2. The left column corresponds 
to models with compactness C = 0.12, while the right column corresponds to C = 0.18. The resolution is Hs2d from Table II. A comparison 
is made with the results presented in [1], where similar solutions were obtained from the spectral code LORENE [112]. 


Type Patch 
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O 

d 
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10^ 
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96 

96 

12 
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o 

d 
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10^ 
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96 

96 

12 
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- 
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- 

32 

- 

40 
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96 

96 
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TABLE III. Grid parameters used for the irrotational solution of compactness C — 0.18, presented in Fig. 7. The separation between the two 
neutron stars is ds/fs = 2.5/0.76 = 3.29. 

























































































































































16 



2 

1 



-1 - : 

-2 .... 

-10 12 3 

X 



X 



X 




FIG. 7. Irrotational binary solution with Vs = 0.76 and compactness C = 0.18. The separation between the two neutron stars is ds/ra = 
2.5/0.76 = 3.29. Left column: contour plots of the conformal factor ip from 1.0 to 1.33 with step 0.01, of the rest-mass density from 0 to 0.3 
with step 0.01, and of the velocity potential 'F from —0.2 to 0.2 with step 0.01, all on the {x, y) plane. Right column: shift and fluid velocity 
vector fields on the {x, y) plane, and contour plot of the rest-mass density on the {y, z) plane. Note that the green sphere corresponds to the 
excised sphere 5'e of COCP-1. 
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ds/rs 


FIG. 8. Top panel: relative error in the relation Mj, = for 

irrotational sequences with C = 0.18 using the Hs2d grid of Table II 
and the Hs3d grid of Table III. Bottom panel: relative error in the re¬ 
lation = Q.dJ for the irrotational sequences of compactness 

C = 0.18 as a function of the coordinate separation ds/rg. 





FIG. 9. Relative change in the central rest-mass density for the irrota¬ 
tional sequences with C = 0.18 of Fig. 6 as a function of separation. 
Note that the decrease as the stars approach each other is 1 order of 
magnitude (or more) smaller than for corotating binaries (cf.. Fig. 3). 


As done in Sec. IV A for corotating binaries, we compare in 
Fig. 6 our irrotational solutions for compactness C = 0.12 and 
0.18 against the corresponding results presented in Ref. [1]. 
Although we use here a relatively small resolution, i.e., Hs2d 
from Table II, the relative difference with Ref. [ 1 ] is again of 
the order of 0.05%. Note that for binaries with compactness 
C = 0.12, the variable Vg ranges from = 0.5 to = 
0.87, which corresponds to coordinate separations dg/rg = 
2.5/0.5 = 5 and dg/vg = 2.5/0.87 = 2.87, respectively. On 
the other hand, for C = 0.18, Vg varies from Vg = 0.5 to 
rg = 0.79, which corresponds to separations from dg/vg = 5 
to dg/vg = 3.16, respectively. Note also that the minimum in 
these plots marks the mass-shedding limit and the creation on 
the equatorial plane of a cusp in the rest-mass density. 

In Fig. 7 we present the results relative to the irrotational 
binary solutions with rg = 0.76 and compactness C = 0.18. 
More specifically, in the left column we report the contour 
plots of the conformal factor ip from 1.0 to 1.33 with step 
0.01, of the rest-mass density from 0.0 to 0.3 with step 0.01, 
and of the velocity potential $ from —0.2 to 0.2 with step 
0.01, all on the (x, y) plane. On the other hand, in the right 
column we show the shift and the fluid velocity vector fields 
on the (x, y) plane, and a contour plot of the rest-mass density 
on the (y, z) plane. 

Similarly, in Fig. 8 we report two global error indica¬ 
tors computed for irrotational binaries with compactness C = 
0.18. More specifically, the top panel shows the fractional 
difference in the Komar and ADM masses, while the bottom 
panel shows the fractional error of the = DdJ rela¬ 

tion; note that the latter is a rather stringent test and that a frac¬ 
tional error below 0.7% gives us confidence on the accuracy 
of our solutions already at an intermediate resolution. Finally, 
in Fig. 9 we plot the relative change in the central rest-mass 
density as the coordinate separation between the two stars is 
reduced. This figure should be compared with the correspond¬ 
ing Fig. 3 for corotating binaries and shows that again the cen¬ 
tral rest-mass density decreases as the two stars approach, but 
also that this decrease is smaller, of 1 order of magnitude or 
more, than in the corotating case. 


C. Spinning sequences 

We conclude our discussion of the results with the new 
COCAL by presenting our first calculations of quasiequilib¬ 
rium binary systems of spinning neutron stars. The neutron- 
star matter is modeled using a piecewise polytrope represen¬ 
tation of the APRl EOS [104]. As mentioned in Sec. IIC, an 
EOS with N polytropic segments requires 2N parameters to 
be specified, which can be thought of as one adiabatic con¬ 
stant, N — \ dividing rest-mass densities, and N adiabatic 
indices. In Ref. [101] it was found that a number of tabu¬ 
lated nuclear matter EOS can be modeled with three segments 
above nuclear density and one in the crust, thus with a total 
of four polytropic zones. The error in the approximation is 
~ 0.1%, or at worst ~ 4%. A fit with a minimum error was 
described in [101] that had a fixed crust with Fq = 1.35692, 
Kq = 3.59389 X 10^^, and three core zones with adiabatic ex- 
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Q. 



FIG. 10. Top: dimensionless binding energy (left panel) and angular momentum (right panel) as a function of the dimensionless orbital 
frequency for sequences of constant rest mass neutron-star binaries with Mo — 1.5388 Mq. Shown with different lines are an: irrotational 
sequence (violet solid line) and two spinning ones with A = 0 and either So = 0.01 (red solid line) or So = 0.05 (blue solid line). All binaries 
are modeled with the APRl EOS and are also shown for comparison is the 4PN irrotational (black dashed line), and the 3PN corotating (green 
dotted line) approximation. Bottom: Fractional differences of the So = 0.01 and So ~ 0.05 spinning sequences with respect to the irrotational 
sequence, AE/Eir := Flsp/i?ir — 1, AJ/Jir := Jap/Jir — 1. 
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14.700 
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log pa 

log Pc 

log Cc 

So 

A 

(Jsp - Jir)/2 

1.666 
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2.294 

7.995 

0.307 
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36.340 

15.642 

- 

- 

- 

1.350 

1.539 

1.636 

9.138 

0.218 

15.221 

35.562 

15.276 

- 

- 

- 

2.661 

1.539 

1.636 

10.246 

- 

15.220 

35.559 

15.275 

0.01 

0 

0.0225 

2.661 

1.539 

1.636 

10.249 

- 

15.220 

35.560 

15.275 

0.01 

-6 

0.0062 

2.661 

1.539 

1.636 

10.260 

- 

15.218 

35.556 

15.274 

0.05 

0 

0.1176 

2.661 

1.539 

1.636 

10.248 

- 

15.219 

35.559 

15.275 

0.05 

-6 

0.0313 


TABLE IV. Top part: summary of the parameters in the piecewise polytropic EOS used to represent the APRl EOS. Bottom part: the first lines 
report the main properties of the maximum-mass nonrotating configuration of the APRl EOS; the second line reports the properties of the 
nonrotating configuration used to construct the constant rest-mass and spinning binary sequences computed in Sec. IV C. The last four lines 
report the stellar properties for the binary in the sequence having the smallest separation, with the last column providing a possible estimate of 
the maximum absolute contribution of the spin angular momentum, (Jsp — Jir)/2. 


ponents {Fi, r 2 , Fa}, joining the different pieces at rest-mass 
densities pi = gr/cm^, and p 2 = 10^® gr/cm^. Ad¬ 

ditional information on the properties of the initial data are 
collected in Table IV. 

The spin contribution to the fluid velocity is expressed 
through the spatial three-vector P [cf., Eq. (56)], which we 
express as 

P = So{-y,x,0), (126) 

where the Cartesian coordinates x, y are centered in the COCP 
patch, and the positive (negative) constant Sq denotes the 
magnitude of corotation (counter-rotation). 

In Table IV we report the properties of a sequence of bi¬ 
nary neutron stars with constant rest mass Mq = 1.5388 Mq, 


corresponding to an ADM mass = 1.35 Mq when the 

stars are at infinite separation. The freedom in the choice of 
the spin velocity vector defined in Eqs. (56) and (126) has 
been fixed by taking A = 0, while for Sq we examine two 
cases: Sq = 0.01, which corresponds to a spinning period 
of ~ 3 ms, and Sq = 0.05, which corresponds to the ex¬ 
treme case of a period ^0.6 ms. These choices correspond to 
spinning periods that are more than twice smaller than those 
considered in [70], where the maximum value considered was 
6.7 ms. Although it is rather unlikely that such small rotation 
periods are encountered in reality in binaries about to merge, 
it is a good consistency check for our new code and an explo¬ 
ration of its limits. 

The top left panel of Eig. 10 reports the dimensionless bind- 














































































































19 



Q. 


FIG. 11. Spin contribution to the total angular momentum. Shown 
for the same binaries presented in Fig. 10 but with A = —6, is the 
relative difference between the angular momentum of the spinning 
binaries and that of an irrotational binary. Note that despite the short 
spin periods considered, the contribution of the spin angular momen¬ 
tum is at most 1% of the total, smaller than the A = 0 case Fig. 10, 
where the same contribution was close to 4%. 


ing energy Eh /Moo of the binary as a function of the dimen¬ 
sionless orbital frequency flMoo- Considered and compared 
are an irrotational binary (violet solid line) and two spinning 
binaries, one with Sq = 0.01 (red solid line) and another 
one with Sq = 0.05 (blue solid line). All binaries are mod¬ 
eled with the APRl EOS using the grid parameters of Table 
V, and both of the spinning binaries have velocity field with 
A = 0. Also shown for comparison is the irrotational fourth 
post-Newtonian (4PN) (black dashed line) as well as the third 
post-Newtonian (3PN) corotating (green dotted line) approx¬ 
imation [125-128]. Explicit forms for these curves are given 
in given in Appendix E. In the top right panel of Eig. 10 we re¬ 
port instead the analogue curves for the dimensionless angular 
momentum J/M^. 

We also note that a closer comparison with a PN expression 
for spinning neutron stars would be interesting. At the same 
time, it involves a number of subtleties. In fact, in order to 
correctly plot such PN spinning curves it is necessary to take 
into account two different effects. Eirst, at each separation 
(orbital frequency) the PN expressions should use the correct 
values of the spins, which we recall increase (in our case ap¬ 
proximately linearly) along the constant rest-mass sequence. 
This is made difficult by the fact that these spins cannot be 
measured separately in our approach as they are part of the 
global solution. These additional terms for aligned spins will 
typically move the binding energy to more negative values 
relative to the irrotational curve. Second, the masses appear¬ 
ing in the PN expression need to modified to account for the 
spin kinetic energy. In contrast to binary black holes, where it 
is possible to distinguish the irreducible mass from the spin- 
induced mass [125-128], accounting for this contribution is 
not easy for neutron stars. However, what is important here is 
that these terms are positive and are will move the quasiequi¬ 
librium curve upwards relative to the irrotational curve and 


toward the corotating solution. 

This is indeed the behavior that is shown by our solutions, 
which fall between the irrotational sequence and the coro¬ 
tating sequence; furthermore, the use of larger initial spins 
yields binding energies that are systematically larger as part of 
the orbital kinetic energy is channeled into “spinning-up” the 
stars. Clearly, the differences in the binding energy with re¬ 
spect to the irrotational binaries are very small even for these 
high spinning rates, with a maximum of 0.03% as it can be 
seen in the bottom left panel of Pig. 10. More pronounced are 
the differences for the angular momentum which are depicted 
on the bottom right panel of Pig. 10 with a maximum of 4%. 
Note also that the sign of 5*0 determines the relative position of 
the spinning sequence relative to the irrotational one. In par¬ 
ticular, with a negative value for Sq, the angular-momentum 
curve for the spinning sequence would have appeared below 
the irrotational one. 

Pinally, in Pig. 11 we report an estimate of the spin con¬ 
tribution to the total angular momentum for a different value 
of A. Por the same binaries presented in Pig. 10 the relative 
difference between the angular momentum of the spinning bi¬ 
naries, Jsp, and that of an irrotational binary, Jq. with A = —6 
is at most 1% of the total, smaller than the 4% value obtained 
with A = 0 of Pig. 10. Note that this quantity is not expected 
to be constant along this sequence, where only the rest mass 
and the spinning coefficient Sq are kept constant. 

If we estimate the dimensionless spin angular momentum 
to be 


X := 


M2 

ADM 


'^sp J'n 
M2 

ADM 


(127) 


then because the ADM mass of the spinning binaries is very 
close to the irrotational one and lies in the range ^ 

[1.33,1.34], the dimensionless spin takes values in the range 
X G [0.027,0.066] for the case with A = 0 and Sq = 0.05, 
and values in the range x G [0.0017, 0.0035] for the A = —6, 
Sq = 0.01 case. In all the cases considered we have found 
that the error estimates discussed in the previous sections lead 
quite generically to relative errors that are < 0.7%, and that 
smaller errors can be obtained with grids having a higher res¬ 
olution. 


V. CONCLUSION 

We have presented the extension of the COCAL code to 
treat binary configurations of compact stars within the IWM 
formalism of general relativity. As with the work done for bi¬ 
nary black holes, we have used multiple coordinate patches so 
as to be able to treat asymmetric binaries. Also in the spirit of 
previous work, we have introduced a particular normalization 
scheme that allows us to accurately compute binary systems 
that have small or large separations, recovering the spherical 
limit at large distances. This is done by keeping the stars at 
fixed coordinate positions, but artificially reducing their ra¬ 
dius. Purthermore, we have made use of surface-fitted coordi¬ 
nates to describe accurately the stellar shape as it varies along 
sequences of constant rest mass. 
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72 

12 
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- 

24 

- 

30 
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TABLE V. Grid parameters used for the irrotational and spinning sequences with the APRl EOS, presented in Eig. 10. 


Also for the nonvacuum spacetimes considered here, we 
have employed the KEH method [76-81], in which the grav¬ 
itational equations are solved using Poisson solvers with ap¬ 
propriate Green’s functions, while for the conservation of rest 
mass we employ a least-squared algorithm. The code makes 
use of a piecewise polytropic description to represent the EOS 
of stellar matter and, for the specific cases considered here, we 
have adopted the representation of the APRl EOS [104]. 

Making use of a suitably adapted formulation described 
in Ref. [63], the code is able to describe fluid flows within 
the stars that corresponds to corotating, irrotational, but also 
spinning binaries. As a validation of the numerical solutions, 
we have constructed a number of sequences of corotating 
and irrotational binary neutron stars having the same mass. 
The results for corotating and irrotational binaries have been 
compared with those published from the pseudospectral code 
LORENE [1], and they revealed that the relative difference in 
the results between the two codes is of the order of 0.05%, 
even when a medium resolution is used for COCAL . 

When considering spinning binaries, and although the code 
can handle arbitrary rotation prescriptions for the individual 
stars, we have concentrated here on the case of fluid flows in 
which the spins are parallel to the orbital angular momentum. 
Eor this class of solutions, and to explore the possible range of 
behaviors, we have considered sequences with stars that either 
are slowly spinning or are spinning at rates that are 10 times 
larger than those observed in binary pulsars systems. In all 
the cases considered, we have found that error estimates of 
different types lead to relative errors that are < 0.7%. 

A number of applications of these results and of additional 
developments of the code are expected to take place in the 
coming months. Eirst, we will explore the impacts of stellar 
spins in numerical simulations of binary neutron stars; more 
specifically, by exploiting the high convergence order of our 
new numerical general-relativistic code [129], we plan to ex¬ 
tend the work carried out in [130] for the inspiral part and the 
one recently published in [30, 34] for the postmerger signal. 
Second, by combining the approaches followed in the solution 
of binary black holes and binary neutron stars, we will extend 
the code to handle also binaries comprising a black hole and 
a neutron star of different masses and spin orientation. Third, 
we will explore the space of solutions in which the spins of the 
neutron stars are oriented arbitrarily as these are likely to cor¬ 
respond to the most realistic configurations. Einally, working 
on a parallelization of the code will allow us to obtain results 
with much smaller computational costs, enabling us to pro¬ 
vide public initial data for spinning binary neutron stars under 
a variety of conditions. 
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Appendix A: Mass and angular momentum 


In this appendix we review the mathematical definitions of 
several of the quantities that have been used to characterize the 
properties of the binaries. We start with the rest mass of each 
star. Mo, defined as an integral over the spacelike hypersur¬ 
face St of the rest-mass density as measured by the comoving 
observers 

Mq := f pu°‘dSa = [ pu°‘yat\/—g(Px 

= / pu^aip^r"^ sin 6drdOdcj). (Al) 

7 St 

In COCAL , integrals like this are computed in dimensionless 
form using normalized coordinates. With the help of Eq. (66), 
Eq. (Al) is rewritten as 

Mq := Rq [ Kl^^^~^'\^/^^'~^">v}'a'tp^f^sin9dfd9d(i), 
7st 

(A2) 

where Ki depends on f. The integrand in Eq. (A2) is evalu¬ 
ated on the gravitational coordinates therefore an interpolation 
to the surface-fitted coordinates is needed before the integral 
evaluation. 

Next, a measure of the total energy of the system is given by 
the ADM mass, which is defined as a surface integral 

at spatial infinity as 

- rr^djlmndS. 

= -^ [ d^'tljdS, = sinedOdcj), 

27r 27r dr 

(A3) 


and which in normalized coordinates becomes 

Madm-=-^[ sin 0d0d(j). (A4) 

J r—Vb 
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Note that spatial infinity in COCAL is represented by a spher¬ 
ical surface with radius r Ri O.Sr^ of the ARCP coordinate 
patch. Closely related to the ADM is the Komar mass of the 
binary, which is related to the timelike Killing field and is 
defined as 


- -4^ 


J = ^j d^dS^ , (A5) 


or in normalized form as 


Mk := 


i?0 

47r 



f ^ sin 9d9d6. 
or 


(A6) 


The angular momentum of the system is also calculated 
from a surface integral at spatial infinity 


J :=-!- / dSi = ^ [ Aij(l>’xWaoSiii9d9d(j), 
Stt Joo Stt 7^ " 


(A7) 


where cj)'' is the generator of the orbital trajectories and we 
have used the maximal slicing gauge. The corresponding nor¬ 
malized quantity is 

-^0 / Aij ^x^fb sin dd9d4>. (A8) 

oTT J 


the interior, instead, we need to solve the TOV equations 


dA 2A dp 
dr e + p dr ^ 
dp (e + p){m + Airr^p) 
dr r^ — 2mr 


(B3) 

(B4) 


where 


dm 

dr 


47rr^e, 


and 


B{r) 


1 

1 — 2m{r)jr 


(B5) 


Of course it not difficult to derive the TOV equations in the 
isotropic coordinates (B2) and then perform a direct numeri¬ 
cal integration in these coordinates. However, this is not nec¬ 
essary and it is possible to always work in Schwarzschild co¬ 
ordinates rescaling the radial profile of the solution so as to 
make the surface of the star appear at the correct position and 
automatically obtain a smooth solution in [0, oo) without re¬ 
sorting to a postprocessing rescaling. 

Comparing Eqs. (Bl) and (B2) it is easy to deduce that 


ip'^{r)dr = \fBdr, and ip‘^{r)r = r , (B6) 


which yield 


Finally, we also compute the “proper mass” of each star as the 
integral of the total energy density measured by the comoving 
observer 



(B7) 


Mp := 


eu'^dSa 


(A9) 


Appendix B: Isotropic coordinates TOV solver 

In this appendix we describe our implementation for ob¬ 
taining spherical solutions and the related rescaling that is 
used in COCAL . We can obtain the same solutions using 
a one-dimensional KEH solver that mimics the full three- 
dimensional code in a 3 h- 1 setting. However, because most 
of the time the TOV equations are presented in terms of 
Schwarzschild coordinates while the actual calculations are 
performed in isotropic coordinates, in what follows we show 
how to transform the system of equations from Schwarzschild 
to isotropic coordinates without having to go through a new 
derivation of equations and automatically obtaining a smooth 
solution at the stellar surface. The results are of course iden¬ 
tical to machine precision, at least for simple polytropes we 
have checked. To the best of our knowledge this approach has 
not been presented before in the literature. 

We recall that the line element in Schwarzschild and in 
isotropic coordinates is given, respectively, by 


Using Eq. (B7), we can rewrite the TOV system in terms of 
the isotropic radial coordinate r as 


dm 

dr 

dp 

dr 

dip 

df 



{e + p){m + Airr^p) 1 
- 2m/r rf ’ 



da m + Airr^p 1 

— = a .- r , 

dr y'l - 2mlr rr 


(B8) 

(B9) 

(BIO) 

(BID 


where we used Eq. (B6) to derive Eqs. (BIO) and (Bll)- B 
is possible to simply integrate the system above to obtain the 
star profile in isotropic coordinates. Initial values at f = 0 are 
needed and although this is not a problem for r, m, and p (the 
latter being in general a free parameter), but the values of a, ip 
are not available to have a smooth matching at the surface of 
the star f = R, whose position is still unknown. On the other 
hand, one way to obtain a smooth solution across the star’s 
surface is to exploit the coordinate transformations 


ds^ = -A{r)dR -b B{r)dr^ + r^dff , 

_ 7-2 , -2 7o2 


(Bl) 

), (B2) 

r = r 


(B12) 

functions 
star. For 

T 

^ ~ 2 

/ 2M M 

1 + y 1 

\ r r 

, (B13) 
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FIG. 12. Solution of isotropic TOV using correct boundary condi¬ 
tions. A rescaling was done at the surface to be at r = 0.625. 


to derive the expressions of ip, a in the Schwarzschild coor¬ 
dinates 


V'(r) = 


r 

M 




(B14) 


a{r) = 



(B15) 


Again making use of Eq. (B7) and of the analytic integration 
of its left-hand side, Eq. (BIO) can be written as 


dint/) = F{r, m{r))dr , 


where 


F(r^ m{r)) 



•\/r^^'2rn^r)7r 


(B16) 


(B17) 


As a result, an integration between r = 0 and r = R of 
F{r,m) leads to the following condition for the conformal 
factor Ip 


StJj” up to ’’Invert V^, compute ip”) have been described in 
detail in [73]. 

The code first initializes the lapse function and the con¬ 
formal factor from some initial spherical solution. Eor 
that purpose we have developed two methods. One is an 
isotropic TOV solver (see Appendix B), and another is a one¬ 
dimensional KEH method. This last choice reproduces the 
KEH approach used in three-dimensional computations, but 
in a one-dimensional mesh. Included in this method are all the 
important ingredients of the three-dimensional code, such as 
the renormalization of variables. Comparing the results from 
these two independent schemes gives us confidence about the 
robustness of the COCAL iterative solutions. 

After a choice of the velocity fluid potential, of the orbital 
angular velocity, and in the case of spinning binaries, also of 
the rotational states of each compact object, the code proceeds 
to the main part of the iteration, which always starts by inter¬ 
polating q = p/p, di^, and s® from the surface-fitted coordi¬ 
nates to the gravitational coordinates. The interpolated quan¬ 
tities are then used in the gravitational Poisson solver, which 
is executed in addition to the root-finding routine explained 
in Sec. IIIB. As discussed there, the constants related to the 
Euler integral, the orbital angular velocity, and the scaling of 
our grids C, Ct, Rq, are calculated at this point, and the lapse 
function, as well as the conformal factor, are updated accord¬ 
ing to 

.^new ^ ^new ^ ^ 

(Cl) 

When the gravitational solver ends, ip,l3‘^,a are interpo¬ 
lated to the surface-fitted coordinates in preparation for the 
fluid Poisson solver. The main steps now are the computation 
of the new value of q, by the use of Eq. (49) or Eq. (32), and 
then the solution of the conservation of rest mass, Eq. (57). 
At this point, also the surface of the star is computed. 

At each iteration step, the fluid computation is performed 
a few times (4 times for the results presented here) since this 
results in a more stable final computation. A relaxation pa¬ 
rameter ^ is used when updating a newly computed variable. 
If is the n-th step value, and $(a;) the result of the 

Poisson solver, then the (n + l)-th step value will be 


ip{0) = ip{R) exp 


F{r, m{r))dr 


(B 18 ) 


so that Ip (and similarly a) at the star’s surface are guaranteed 
to match smoothly the exterior solution, as it can be seen in 
Pig. 12. 


Appendix C: Iteration scheme 

The iteration procedure for binary stars is similar to the one 
for binary black holes described in Ref. [73], but it also con¬ 
tains the fluid coordinates, and a solution for the extra fluid 
variables in a multipatch setting. An overall picture of this 
procedure is shown in Pig. 13, where the steps of the gravita¬ 
tional Poisson solver (essentially everything from ’’Compute 


$(®®+i)(a;) .= + (1 _ ^)$(")(a;), (C2) 

where 0.1 < ^ < 0.4. Usually ^ = 0.4 for a,ip,P'^,q, while 
^ = 0.1 for p. 

The criterion used by COCAL to stop the iteration is given 
by 

for all points of the grids, and all variables a, ip, /3®, q, 4>, 
where we used ec = 10“® in this paper. In almost all of our 
calculations, ip and a converge to machine precision, while 
the error in the fluid variables q and $ decreases to « 10“^^ 
before the error in the shift reaches 10“^. This is due to the 
existence of points in the gravitational mesh where the shift 
has almost zero values, and convergence is much slower there. 
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FIG. 13. The COCAL iteration for spinning binary stars. 
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Type Patch 

Ta 

Ts 

rg 

Tc 

re 

n!^ 

A/ 

A” 

Nr 

Ng 

N^ 

L 

WD COCP 

- 1 

0.0 

1.0 

10^ 

1.50 

1.25 

64 

64 

96 

192 

144 

144 

12 

COCP 

- 2 

0.0 

1.0 

10^ 

1.50 

1.25 

64 

64 

96 

192 

144 

144 

12 

ARCP 


5.0 

- 

10® 

6.25 

- 

16 

- 

20 

192 

144 

144 

12 


TABLE VI. Grid parameters used for the white-dwarf solutions with T = 5/3. 



FIG. 14. Binary white dwarfs solution with compactness C = 2 x 10“^, stellar centers placed at a: = ±1.25 and unit radii. Shown from 
left to right are: a contour plot of the lapse function from 0.9994 to 1.0 with step of 2 x 10~®, the shift vector field with maximum value 
8.4 X 10~®, and a contour plot of the rest-mass density from 2 x 10“® to 10“"^ with step 8 x 10“®. Note that the green sphere corresponds 
to the excised sphere Se of COCP-1. 


Neglecting such points can speed up a solution by a factor of 
at least 2. Currently, COCAL is running on a serial processor 
and it needs around 3-4GB of RAM to produce the solutions 
presented in this work. With an Intel Xeon 3.60 GHz proces¬ 
sor, about two days are needed for these computations, with 
the irrotational configurations taking longer than the spinning 
ones; this is not a surprise since convergence is faster for coro¬ 
tating binaries. 

Appendix D: Corotating binary white dwarfs 

To test the sensitivity of our code and to prepare for future 
work concerning neutron star-white dwarf or black hole-white 
dwarf binaries, we also compute a corotating binary white 
dwarf solution. Here the fields are orders of magnitude less 
than the ones encountered in typical binary neutron star bi¬ 
naries and greater resolution is required in order to acquire 
smooth solutions. The resolution used is reported in Table VI 
where we can see that an increase in Ng, by a factor of 3 
relative to the solutions obtained in Fig. 2, Table II, was used. 
In Fig. 14 we show a representative binary white dwarfs so¬ 
lution with compactness C = 2 x 10“^, with centers placed 


at a; = ±1.25 and unit radii. From left to right, the differ¬ 
ent panels report the contour plot of the lapse function from 
0.9994 to 1.0 with step of 2 x 10“®, the shift vector field, 
and the contour plot of the rest mass density. Note that the 
plots are centered at the origin of the COCP-1 patch and the 
green circle refers to the excised sphere Se- Values inside Se 
are taken from the COCP-2 patch. The shift vector in binary 
white dwarfs is approximately 4 orders of magnitude smaller 
than the one typically encountered in neutron stars, while the 
quantity | a — 11 is about 3 orders of magnitude smaller. Over¬ 
all we see a good convergence between the different coordi¬ 
nate systems even for these small values of the metric quanti¬ 
ties. 


Appendix E: Post-Newtonian approximation 

The 4PN approximation for the binding energy of a system 
of two nonspinning bodies with masses Mi, M 2 and in qua¬ 
sicircular orbit has been used in Fig. 10 to compare with the 
numerical results of irrotational and spinning binaries. The 
explicit expression is given by [125-127] 
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where 7 ^ is the Euler constant, M := Mi + M 2 is the total 
mass of the system, v := MiM 2 /M'^ is the symmetric mass 
ratio and x is the dimensionless orbital frequency 




(E2) 


Eor a corotating binary the binding energy also includes terms 
from the kinetic energy of the spins as well as from the spin 
orbit interaction. To 3PN this extra contribution is 


Ait^coi 

Mc 2 




= (2 — 6z/)x'’ + (—lOi^ + 2bu‘^)x 


(E3) 


and therefore the total binding energy is given by the sum of 
(El) and (E3). 
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Eor a system that satisfies the first law of binary mechanics 
and has binding energy of the form 


N 


Eb{x) = + Bi\nx)x'/‘^ 

i=l 

the angular momentum will be 

N 

Ji.x) = ^ 


(E4) 


2=1 ^ 


* -XAi + Biinx)- 


X 2 


(E5) 


Using (E5) the irrotational part of the angular momentum is 
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while the corotating additional part is 

= (4 - Ylv)x^l'^ + (-16i2 + 40 V)x'^/2 . 
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